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Abstract: In this work, the dependence of reversal potentials and zero-current fluxes on diffusion 
coefficients are examined for ionic flows through membrane channels. The study is conducted for 
the setup of a simple structure defined by the profile of permanent charges with two mobile ion 
species, one positively charged (cation) and one negatively charged (anion). Numerical observations 
are obtained from analytical results established using geometric singular perturbation analysis of 
classical Poisson-Nernst-Planck models. For 1:1 ionic mixtures with arbitrary diffusion constants, 
Mofidi and Liu [arXiv:1909.01192 ] conducted a rigorous mathematical analysis and derived an 
equation for reversal potentials. We summarize and extend these results with numerical observations 
for biological relevant situations. The numerical investigations on profiles of the electrochemical 
potentials, ion concentrations, and electrical potential across ion channels are also presented for the 
zero-current case. Moreover, the dependences of current and fluxes on voltages and permanent 
charges are investigated. In the opinion of the authors, many results in the paper are not intuitive, 
and it is difficult, if not impossible, to reveal all cases without investigations of this type. 

Keywords: Reversal potential, effects of diffusion coefficients, permanent charge 



1. Introduction. 

Ion channels are proteins found in cell membranes that create openings in the membrane to allow 
cells to communicate with each other and with the outside to transform signals and to conduct tasks 
together [3,12]. They have an aqueous pore, that becomes accessible to ions after a change in the 
protein structure that makes ion channels open. Ion channels permit the selective passage of charged 
ions formed from dissolved salts, including sodium, potassium, calcium, and chloride ions that carry 
electrical current in and out of the cell. 

To unravel how ion channels operate, one needs to understand the physical structure of ion 
channels, which is defined by the channel shape and the permanent charge. The shape of a typical 
ion channel is often approximated as a cylindrical-like domain with non-uniform cross-sectional area. 
Within a large class of ion channels, amino acid side chains are distributed mainly over a "short” and 
"narrow" portion of the channel, with acidic side chains contributing permanent negative charges and 
basic side chains contributing permanent positive charges. 

The spatial distribution of side chains in a specific channel defines the permanent charge of 
the channel. The spatial distribution of permanent charge forms (most of) the electrical structure of 
the channel protein. The spatial distribution of mass forms the structure studied so successfully by 
molecular and structural biologists. Ions that move through channels are often only an angrom or 
so distant from the permanent charges residing on acid and base side chains. And electrical forces 
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are in general much stronger than entropic forces. So in most cases the electrical structure is more 
important in determining how ions go through a channel than the mass structure. Sometimes the 
dielectric properties ('polarization') of the protein contribute charge that is significant. Then, the spatial 
distribution of dielectric properties becomes an important part of the electrical structure. 

The most basic function of ion channels is to regulate the permeability of membranes for a given 
species of ions and to select the types of ions and to facilitate and modulate the diffusion of ions across 
cell membranes. At present, these permeation and selectivity properties of ion channels are usually 
determined from the current-voltage (I-V) relations measured experimentally ([12,19]). Individual 
fluxes carry more information than the current, but it is expensive and challenging to measure them 
([23,27]). The I-V relation defines the function of the channel structure, namely the ionic transport 
through the channel. That transport depends on driving forces expressed mathematically as boundary 
conditions. The multi-scale feature of the problem with multiple physical parameters allows the 
system to have great flexibility and to exhibit vibrant phenomena/behaviors - a great advantage 
of "natural devices" ([11]). On the other hand, the same multi-scale feature with multiple physical 
parameters presents an extremely challenging task for anyone to extract meaningful information from 
experimental data, also given the fact that the internal dynamics cannot be measured with present 
techniques. The general inverse problem is challenging, although specific inverse problems have been 
successfully solved with surprisingly little difficulty using standard methods and software packages 

(M). 

To understand the importance of the relation of current and permanent charges, that is, the I-Q 
relation, we point out that the role of permanent charges in ionic channels is similar to the role of 
doping profiles in semiconductor devices. Semiconductor devices are similar to ionic channels in 
the way that they both use atomic-scale structures to control macroscopic flows from one reservoir 
to another. Ions move much as quasi-particles move in semiconductors. In a crude sense, holes 
and electrons are the cations and anions of semiconductors. Semiconductor technology controls the 
migration and diffusion of quasi-particles of charge in transistors and integrated circuits. Doping is 
the process of adding impurities into intrinsic semiconductors to modulate its electrical, optical, and 
structural properties ([41,47]). In a crude sense, doping providers the charges that acid and base side 
chains provide in a protein channel. 

Ion channels are almost always passive and do not require a source of chemical energy (e.g., ATP 
hydrolysis) in order to operate. Rather they allow ions to flow passively driven by a combination of the 
transmembrane electrical potential and the ion concentration gradient across the membrane. For fixed 
other physical quantities, the total current X — I(V,Q) depends on the transmembrane potential V 
and the permanent charge Q. For fixed Q, a reversal potential V = V r ev(Q) is a transmembrane potential 
that produces zero current X( V r ev ( Q), Q) =0. Similarly, for fixed transmembrane potential V, a reversal 
permanent charge Q — Q rei ,(V) is a permanent charge that produces zero current X(V, Q rep (V)) = 0. 

The Goldman-Hodgkin-Katz (GHK) equation for reversal potentials involving multiple ion species 
([20,22]) is used to determine the reversal potential across ion channels. The GHK equation is an 
extension of the Nernst equation - the latter is for one ion species. The classical derivations were 
based on the incorrect assumption that the electric potential <1>(X) is linear in X - the coordinate along 
the longitude of the channel. This assumption is particularly unfortunate because it is the change in 
the shape of the electrical potential <1>(X) that is responsible for so many of the fascinating behaviors 
of transistors or ionic systems ([9,10,37,42,44,46]). There was no substitute for GHK equations until 
authors of [15,38] recently offered equations derived from self-consistent Poisson-Nernst-Planck (PNP) 
systems, to the best of our knowledge. 

There have been great achievements in analyzing the PNP models for ionic flows through ion 
channels ([13,27-29,34,35], etc.). Although mathematical analysis plays a powerful and unique role to 
explain mechanisms of observed biological phenomena and to discover new phenomena, numerical 
simulations are needed to fit actual experimental data and study cases where analytical solutions do 
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not exist. Furthermore, numerical observations may give clues for more theoretical investigations. 
Indeed, numerical and analytical studies are linked; any progress in one catalyzes work in the other. 

The numerical results, throughout the paper, are gained from the algebraic systems (1), (7), 
(8) and (3), obtained from reduced matching systems of analytical results in [38] and [13]. The 
nonlinear algebraic systems are then solved by the Matlab® function /solve that uses the trust-region 
dogleg algorithm. The trust-region algorithm is a subspace trust-region method and is based on the 
interior-reflective Newton method described in [ 8 ]. Our numerical results indicate that current-voltage 
and current-permanent charge and even zero-current relations depend on a rich interplay of boundary 
conditions and the channel geometry arising from the mathematical properties analyzed in ([13,29,38, 
48]). 

This paper is organized as follows. The classical PNP model for ionic flows is recalled in Section 
1.1 to prepare the stage for investigations in later sections. In Section 2, we study zero current problems 
to investigate the corresponding fluxes, and reversal potentials V r ev • In particular, we compare a special 
case of the reversal potential with the GHK equation. Some other numerical observations are also 
provided to study profiles of relevant physical quantities in Section 2.4. In Section 3, we first recall the 
analytical results in [13] when diffusion constants are also involved. Then numerical observations are 
provided to examine behaviors of current, voltage, and permanent charge with respect to each other in 
some general cases. Some concluding remarks are provided in Section 4. 

1.1. Poisson-Nernst-Planck models for ionic flows. 

The PNP system of equations has been analyzed mathematically to some extent, but the equations 
have been simulated and computed to a much larger extent ([1,6,7,24,26]). One can see from these 
simulations that macroscopic reservoirs must be included in the mathematical formulation to describe 
the actual behavior of channels ([18,40]). For an ionic mixture of n ion species, PNP type model is, for 
k — 1,2 


Poisson: 

Nernst-Planck: 


V ■ (e r (X)e 0 V<I>) - -e 0 (l> s C s + g(X)), 

S —1 

dfC/t + V ■ Jit = 0, —Jk — i~^T> k (X)C k 'V Pk> 


( 1 ) 


where X £ Cl with Q being a three-dimensional cylindrical-like domain representing the channel 
of length L (run — Lx 10~ 9 m), Q(X) is the permanent charge density of the channel (with unit 
1M — lmol/L = 10 3 mol/m 3 ), e r (X) is the relative dielectric coefficient (with unit 1), eo ~ 8.854 x 
10 12 F m 1 is the vacuum permittivity, eg « 1.602 x 10 19 C (coulomb) is the elementary charge, kg rts 
1.381 x 10 23 JK ' 1 is the Boltzmann constant, T is the absolute temperature (T « 273.16 K =kelvin, for 
water); <T> is the electric potential (with the unit V — Volt — JC '), and, for the k-ih ion species, Q, is 
the concentration (with unit M), z k is the valence (the number of charges per particle with unit 1), u k 
is the electrochemical potential (with unit J = CV) depending on <E> and Q.. The flux density J k (X) 
(with unit mol/s) is the number of particles across each cross-section in per unit time, T> k (X) is the 
diffusion coefficient (with unit m 2 /s), and n is the number of distinct types of ion species (with unit 1). 

Ion channels have narrow cross-sections relative to their lengths. Therefore, three-dimensional 
PNP type models can be reduced to quasi-one-dimensional models. The authors of [39] first 
offered a reduced form, and for a particular case, the reduction is precisely verified in [33]. The 
quasi-one-dimensional steady-state PNP type is, for k — 1,2,..., n, 


1 d 
MX) dX 


e r (X)e 0 A(X) 


dO 

dX 


= -eo( £> s C s + S(X)) , 


\S=1 J 

s = °- =St d ‘ (x)/1 < x)c ‘®' 


(2) 
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where A(X) is the area of cross-section of the channel over location X. Equipped with system (2), we 
impose the following boundary conditions, for k — 1,2, • • • , n, 

0(0)-V, Cjfc(O) = Lit > 0; O(L) — 0, C k (L) = R k > 0. (3) 


One often uses the electroneutrality conditions on the boundary concentrations because the solutions 
are made from electroneutral solid salts, 

n n 

Y^z s L s — z s R s — 0- (4) 

S—1 S=1 

The electrochemical potential }i k {X) for the k- th ion species consists of the ideal component ji‘ k (X) 
and the excess component jt e k x (X), i.e., }i k (X) — }i l k (X) + ji e k [X). The excess electrochemical potential 
}i e k (X) accounts for the finite size effect of ions. It is needed whenever concentrations exceed say 50mM, 
as they almost always do in technological and biological situations and often reach concentrations 
1M or more. The classical PNP model only deals with the ideal component ji‘ k (X), which reflects the 
collision between ion particles and water molecules and ignores the size of ions; that is, 

MX) = ? 4 rf (X) =z k e 0 <i>(X)+k B Tln Ck ^, (5) 

ho 

where Co is a characteristic concentration of the problems, and one may consider, 

Co = max { L k/ R k/ sup |Q(X)|}. (6) 

l - k - n Xs[0,L] 

For given V, Q(X), L k s and R k 's, if (<J>(X), C k (X), J k ) is a solution of the boundary value problem 
(BVP) of (2) and (3), then the electric current I is 

n 

1 = z s J s . (7) 

S=1 

For an analysis of the boundary value problem (2) and (3), we work on a dimensionless form. Set 


Let 


Vq — max { sup V k (X)} and e r = sup e r (X). 

l<k< n Xe[0,L] 


c 2 _ IreoksT 6(x) _*rW x _X _ A(X) T^X) 

elV Co' r( Er ' V K ) L2 ' Dk[X) - V 0 ' 


0^ 

Q(x) 


QW = ^, *(*)-&!&. 


In terms of the new variables, the BVP of (2) and (3) become, for k — 1,2, ■ • • , n, 


( 8 ) 


= 0 , -J k ^Hx)D k (x)c k ^ x fi k , 


(9) 


with the boundary conditions 


<p{ Q) = V=^=V, Cfc(0) — l k — <p(l) — 0, c k (l)=r k =^. 

k b i Cq Cq 


(10) 
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Remark 1. The actual dimensional forms of quantities have been used for all figures throughout the 
paper, i.e., for k — 1,2, 

C* =C 0 c k (M), Q = C 0 Q (M), y k = (e 0 <i> + k B T\n(C k /C 0 )) (J), 

knT - „ (11) 

<J>=^<£(V), J k = LC 0 V 0 J k (mol/s), 1 = LC 0 V 0 e 0 I (A), 
eo 

and we take Co = 10M, L = 2.5nm and 'Dq — 2.032 x 10 - 9 m 2 /s, and, for diffusion constants ([31]), 
for k — 1,2, 

T> k =1.334 x 10~ 9 m 2 /s for Na + , or 

V k =2.032 xl0” 9 m 2 /s for Cl“, or (12) 

V k =0.792 x 10~ 9 m 2 /s for Ca 2+ . 


1.2. Setup of the problem. 

We give the precise setup of the problem considered in this work. More precisely, we assume 

(A0) The ionic mixture consists oftzvo ion species with valences Zj = —zi = 1; 

(Al) D k (x) = Dj - for k = 1,2 is a constant and e(x) = 1; 

(A2) Electroneutrality boundary conditions (4) hold; 

(A3) The permanent charge Q is piecewise constant with one nonzero region; that is, for a partition 0 < a < 

b <1 of [0,1], 


QW 


Qi = £>3 = 0, X e (0, a) U (b, 1), 
0 , 2 , X e {a, b), 


(13) 


where Q 2 is a constant. 

We assume that e > 0 in system (9) is small. The assumption is reasonable since, if L = 2.5nm = 
2.5 x 10 “ 9 m and Co — 10M, then e « 10~ 3 ([14]). The assumption that e is small enables one to treat 
the BVP (9) and (10) of the dimensionless problem as a singularly perturbed problem. 

A geometric singular perturbation (GSP) framework for analyzing BVP of cPNP systems was 
developed first in [13,34] for ionic mixtures with two types of ion species. This GSP framework was 
extended to arbitrary number of types of ion species successfully only when two special structures of 
the PNP system were revealed ([35]). One special structure is a complete set of integrals (or conserved 
quantities) for the e = 0 limit fast (or inner) system that allows a detailed analysis of singular layer 
component of the full problem. It should pointed out that most of the integrals are NOT conserved 
for the physical problem since, no matter how small it is, £ is NOT zero. The GSP allows one to make 
conclusion about the BVP for £ > 0 small from information of £ = 0 limit systems. The other special 
structure is that a state-dependent scaling of the independent variable turns the nonlinear limit slow (or 
outer) system to a linear system with constant coefficients. The coefficients do depend on unknown 
fluxes to be determined as a part of the whole problem, and this is the mathematical reason for the 
rich dynamics of the problem. As a consequence of the framework, the existence, multiplicity, and 
spatial profiles of the singular orbits - zeroth order in £ approximations of the BVP - are reduced to a 
system of nonlinear algebraic equations that involves all relevant quantities altogether - the physical 
evidence why this framework works since "Everything interacts with everything else" in determining 
ion channel functions. This geometric framework with its extensions to include ion size effects to 
some extent in [28,32,45] has produced a number of results that are central to ion channel properties 
([27,29,30,36,38,48]); for example, it was shown in [29] that a positive permanent charge may enhance 
anion flux as well as cation flux; and, in order to optimize effects of the permanent charge, the channel should 
have a short and narrow neck within which the permanent charge is confined; and, it was shown in [48] 
that, large permanent charge is responsible for the declining phenomenon - decreasing flux with increasing 
transmembrane electrochemical potential. We refer the readers to the aforementioned papers for more 
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details on geometric singular perturbation framework for PNP as well as concrete applications to ion 
channel problems. 

In this paper, we will apply some results and follow the notations in [13] and [38] for analytical 
results where the quantities are all in their dimensionless forms. Besides, for simplicity, we use the 
letters /, r and Q (l where h — h — L r i — ^2 — r , Q = 2 Qq. Also, for definiteness, we choose a — 1/3, 
b — 2/3 and h(x) — 1. 

Remark 2. We remind the readers that the quantities V , /, r, cp, Q, ([>, ftp, Jp, Dp and I are dimensionless 
quantities corresponding to the dimensional quantities V, L, R, Cp, Q, d>, jtp, Jp, Vp and X, respectively, 
obtained from (8). We switch from dimensional form to the dimensionless form and vice versa several 
times throughout the paper. 

The highlights of our studies in this paper as well as in [13,29,38,48] applied to the setup of this 
paper include 

(i) a system for zero-current condition (see system (1)) that can be used to determine the reversal 
potential in terms of other parameters (see (7)); 

(ii) an examination on how the reversal potential depends on permanent charge: its sign and its 
monotonicity in permanent charge (see Section 2.2); and a comparison between this reversal 
potential and that from GHK in the special setting (see Section 2.3); 

(iii) a characterization on monotone dependences of the reversal potential on the ratio of diffusion 
coefficients in terms of different conditions on the boundary concentrations (see Section 2.2), as 
well as effects of un-equal diffusion coefficients on signs of zero-current flux and its dependence 
on permanent charge (see Section 2.1); 

(iv) numerical spatial profiles under the zero current condition of the concentrations and electric 
potential, and hence, the profiles of the electrochemical potentials for several choices of permanent 
charges that reveal special features of permanent charge effects (see Section 2.4, particularly. 
Remark 4); 

(v) numerical and analytical studies of I-V and I-Q relations, and zero-voltage current and its rich 
dependence on permanent charge (see Section 3.3). 

2. Zero current problems with general diffusion constants. 

In this section, we study how boundary concentrations, electric potential, permanent charges and 
diffusion constants work together to produce current reversal. Throughout this section, in order to 
express the effects of diffusion constants on zero-current flux and reversal potential, we study and 
compare the results for different cases of diffusion constants where T>\ — T >2 and where T>\ T> 2 , to 
indicate and emphasize the differences. 

Diffusion is the phenomenon through which the spatial distribution of solute particles varies 
as a result of their potential energy. It is a spontaneous process that acts to eliminate differences in 
concentration and eventually leads a given mixture to a state of uniform composition. Fick's first 
law [16] describes diffusion of uncharged particles by 3fC = T>d \ x c, where c is the concentration, T> 
is the diffusion constant and f is time. Frequently, the determination of diffusion constants involves 
measuring sets of simultaneous values of f, c, and x. These measured values are then applied 
to a solution of Fick's law to get the diffusion constants. Many techniques are available for the 
determination of diffusion constants of ions (charge particles) in aqueous solutions ([2,5,17,31,43], etc.). 
When diffusion constants are equal, classical electrochemistry tells that many electrical phenomena 
"disappear" altogether, e.g., the "liquid junction" is zero. If the diffusion constants of potassium and 
chloride are equal, classical electrochemistry says that KC1 acts nearly as an uncharged species. Indeed, 
that is the basis for the saturated KC1 salt bridge used in a broad range of electrochemical experiments 
for many years. Therefore, the equal diffusion constants case is quite degenerate. Experimental 
measurements are exclusively performed under isothermal conditions to avoid deviation of T> values. 
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Nevertheless, even diffusion constants of certain ionic species may differ from one method to another, 
even when all other parameters are held constant. Many things becomes much more complicated 
mathematically when the diffusion constants are not equal, however. This complexity is what makes 
many biological and technological devices interesting and useful, even valuable. Some kinds of 
selectivity depend on the non-equality of diffusion constants, as well. 

Applying GSP theory to cPNP equations (1) for two ion species with diffusion constants T> k , k = 
1,2, the authors of [38] obtained an algebraic matching system with eleven equations and eleven 
unknowns for zero current problems and singular orbits on [0,1]. They further reduced the matching 
system for the case where two ion valences satisfy zi — — Z 2 - It follows that the reduced matching 
system for zero current I = Ji — J 2 = 0 when zi = —Z 2 = 1 is, 

G\(A, Qo,6) — V and G 2 (A,Q o ,0) = 0, (1) 

where 

G ^- 0) = e (‘" M + A - (I+e >‘" 1 tbf. (2) 

G 2 (A, Qo, 9) —9Q 0 In S s a b l 6 9 Q° o ~ N; 

and, 6 — . Also, A is the geometric mean of concentrations at x — a, that is, 

A = y Jc 1 (a)c 2 {a), (3) 

and 

B — ~~(l ~ A) + r, S a = ^Ql + A 2 , S b = ^Q 2 0 + B 2 , N = A - l + S a - S b . (4) 

In what follows, numerical simulations are conducted with the help of analysis on system (1). The 
combination of numerics and analysis gives a better understanding of the zero-current problems and 
completes the analytical results obtained in [38]. 

2.1. Zero-current flux. 

We aim to clarify the relationships of ion fluxes with permanent charge and diffusion constants 
when current is zero. Recall that fluxes }\ and J 2 are equal for this case and let / denote it. For 
any permanent charge Q — 2Qq, once a solution ( A , V) of (1) is obtained, it follows from matching 
equations (see Appendix in [38]) that / is given by 

6D 1 D 2 (A-l) _ 6DiD 2 (r-B) 

1 (Dx +D 2 ) (D x +D 2 ) ' U 


Sign of zero-current flux J. It was observed in [15] that the Nemst-Planck equation in (9) (with 
dimensionless quantities) gives, for k — 1,2, 


A 

D k 


f 1 1 1 

/ ,— r^dx — z k V + In-. 

Jo h(x)c k (x) r 


( 6 ) 


Therefore, the sign of flux J k depends only on the boundary conditions /, r and V. Note that (6) holds 
for any condition, not just zero-current condition. 

For zero-current problem, V — V rev depends on /, r, D |, D 2 , and Q as well, in general. So the sign 
of zero-current flux / seems to depend on all quantities and to be difficult to figure out. It is not the 
case. A consequence of (5) together with Theorem 3.4 in [38] is that: 
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The zero-current flux J has the same sign as that of l — r. 

The latter follows directly from a statement in Theorem 3.4 in [38] that, for zero-current, / — A 
has the same sign as that of / — r. That is consistent with observations in Figure 1 where T>\ — 
1.334 x 10 9 m * (i) 2 * * * * * * /s is fixed, and 'Do varies from the same value (i.e., p — 1), to Do — 2.032 x 10” 9 m 2 /s 
(i.e., jo > 1), to a random large value (i.e., p 1). 




Figure 1. The function J — J (Q) for various values of p = Do /Dp The left panel for L = 2mM and 

R = 5mM; the right panel for L = 5 mM and R = 2 mM. 

Dependence of zero-current flux J on Q and Dp's. Concerning the dependence of the zero-current 

flux J on Q, we have the following. 

(i) IfD\ — Do, then the zero-current flux J is an even function in Q, and it is monotonically increasing 
for Q > 0. 

In this case, 6 — 0 and hence, it follows from (l)-(3) that A can be determined from N — 0. Since 
N is an even function in Qq, A is an even function in Qq, and so is the zero-current flux J from (5). 

(ii) IfDDo, then the zero-current flux J is not an even function in Q and the monotonicity of the 
zero-current flux J in Q seems to be more complicated. 

In this case, it can be seen that Go in (2) is not an even function in Q, and hence, the zero-current 
flux J is not. We would like to point out that, it follows from [48], for fixed p — Do/D\, no matter 
how large, one always has J -X 0 as Q —? ±oo, that is consistent with the observations in Figure I. 

(iii) Another fascinating result is that the magnitude of p — Do/D\ affects the monotonicity of the 
zero-current flux J in Q. 

In this case, if one fixes D\, and let Do increases from small values to Do —> oo, (i.e., p —> oo), then 
it follows from (5) that there is a meaningful change in the monotonicity of flux, for small values of Q, 
that is not intuitive. Let us consider the case where L < R and Q < 0 is small. It follows from (1) and 

(2) that when Q increases, for the geometric mean of concentrations A, one obtains 

(a) A increases if p ~ 1 (that is 0 « 0), and consequently the zero-current flux J decreases; 

(b) A decreases if p 1 (that is 0 A> 1), and hence, the zero-current flux J increases. 

Thus, depending on the size of p, the zero-current flux J may increase or decrease for Q < 0 

small, which is also consistent with the observations in Figure 1. The analysis for the case with L > R 
is similar. 

It seems likely that the engineer, like evolution, will use these mathematical properties to control 
the qualitative properties of channels, technological and biological. 

2.2. Reversal potential V rev . 

Experimentalists have long identified reversal potential as an essential characteristic of ion 
channels [21,25]. Reversal potential is the potential at which the current reverses direction, i.e., 
V = <b(0) - ^(L) that produces zero current I. Using dimensionless form of quantities (see Remark 
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2), it follows from (1) and (2) (where there are two ion species with valences z\ — —z 2 = 1) that for 
general permanent charge Q = 2Qq f 0 with arbitrary diffusion constants ([38]), the variable A can be 
solved uniquely from G 2 — 0 in (1), and the reversal potential is then 




= fl(ln 


Sfl + 0Qo 

Sk + 0Qo 



(1 + 0) In 


A(Q 0 ,9) 

B(Qo, 0) 


+ In 


S a 


s b 


Qo 

Qo 


(7) 


Range of Reversal potential V rev . For fixed /, r, and for any given Q — 2Qq, it follows from Theorem 
4.2 in [38] that there exists a unique reversal potential V rev such that V rev < | In ( |. As Q —> ± 00 , then 
V rev gets close to the boundary values, i.e. V rev —f ± In (. 

Zero reversal potential. One particular case is when the reversal potential is zero. To examine under 
what conditions one obtains V rev — 0, it follows from Theorem 4.2. in [38] that, 

(i) ifDi = D 2 , then V rev (Q ) =0forQ = 0, 

(ii) ifD] < Do, then there is a Q < 0, such that V rev (Q) — 0, 

(hi) ifD 1 > Do, then there is a Q > 0, such that V rev (Q) — 0. 

Considering the second case above, the observations in Figure 3 show that as p — D 2 /D\ increases, 
magnitude of the corresponding Q becomes larger. In fact, as p -+ 00 , then Q —> — 00 . 

Reversal potential V rev (Q) for Q = 0. For Q ~ 0, one has V rev (0) — Bln 1 - from Theorem 4.2 in [38] 
where 9 — (D 2 — Di)/ (D 2 + Dj). Therefore, 

(i) ifDi — D 2 , then V rev (0) — 0, 

(ii) ifD\ 7 ^ D 2 , then V rev (0 ) has the same sign as that of 9(1 — r). 

Let us consider the case where D 1 < Do. In that case, V rcv ( 0) has the same sign as that of l — r. 
This is reasonable, since for V = 0 we have \Ji\ < \} 2 \ (since all but Jp/ D/ f are independent of Dp in 
(6)), and to help \Ji \ more than I/ 2 1 to get /) = J 2 for zero current condition, one needs to increase V 
when l > r (and decrease V when 1 < r), and that is why V rev (0) > 0 for l > r (and V rev (0) < 0 for 
l < r). This is consistent with observations in Figure 3 as well. The analysis for the other case with 
Di > D 2 is similar. 

Monotonicity of Vrev with respect to Q. It follows from Theorem 4.4 in [38] that, for 0Q > 0, d q V r ev 
has the same sign as that of l — r. This analytical result does not allow conclusions about the case for 
8Q < 0, though. The observations in Figures 2 and 3 show that result holds for any 6 and Q. Thus we 
have 

Conjecture: V rev is increasing in Qfor l > r and decreasing in Qfor l < r. 

We remark that, in Figure 2, we take L — 20mM, R = 50mM, and D\ — 1.334 x H+ 9 rrr/s and 
T >2 — 2.032 x 10 9 m 2 /s which are diffusion constants of, say Na + and Cl - respectively (see the solid 
line), and T>\ — 1.334 x 10~ 9 m 2 /s and T> 2 — 0.792 x 10~ 9 m 2 /s, where T> 2 is the diffusion constants of 
Ca 2+ (see the dashed line). 
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Figure 2. V = V re i>{Q) decreases when L < R, independent of values of diffusion constants. 




Figure 3. The function V = V re v(Q )■' The left panel for L = 2mM and R = 5 mM; the right panel for 

L = 5 mM and R = 2 mM. 

Dependence of V r ev on p — Do/ D\. In terms of the dimensionless form of quantities (See Remark 
2), let us examine the dependence of V rev on p — Di/D\ for effects of D\ and D 2 . It follows from 
Proposition 4.6 in [38] that 

The reversal potential V rev is increasing in p if 1 > r and is decreasing in p if l < r. 

This feature reveals a fantastic aspect that is not intuitive immediately. Recall the equation (6). 
Given the boundary values and diffusion constants, the values one obtains for all terms in (6) except Jp 
are independent of D/ f ([35]). The relation surely holds for the zero-current condition, i.e., /1 = jo with 
V = V rev . Now, let fix D 1 and increase D 2 (so p is increasing). Then I/ 2 1 increases since all but Jo/Do 
in (6) are independent of D 2 . Consequently, to meet zero-current condition, we need to increase \Ji |. 
Intuitively increasing V rev seems to lead to an increase in | J\ |. This intuition agrees with the results 
only for / > r. In the other case, i.e., I < r, it is the exact opposite, though. That is, for l < r, it says, as 
p increases, V rev decreases. This counterintuitive behavior could be explained by the fact that c\{x) 
depends on V rev , and reducing V rev could still increase /) |. In fact, l < r will result in reducing V rea , 
but C\ (x) changes in a way that consequently increases | J\ \. 

To illustrate the counterintuitive behavior, we provide a numerical result in Figure 4. We choose 
Co, L and T>\ for Na + as mentioned in Remark 1. Now, suppose that T>\ — 0.792 x 10 9 m 2 /s, 
and consider the boundary concentrations L — 20 mM, R — 50 mM and Q — 1 M. In this case 
Vrev — —16.7657 mV and J — —10.6184 Mmol/ (s N/\), where Mmol = 10 6 mol. Now, if we increase 
T >2 to T >I = 2.032 x 10“ 9 m 2 /s which is Cl diffusion constant, then V rev — —19.5527mV and J — 
—11.3146 Mmol/ (s N^). These values make sense now, based on the above discussion. Note that, we 
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283 just pictured the middle part of the channel in Figure 4 since the sides are almost identical. One should 

284 notice that it is hard to realize, from Figure 4, how L < R will result in reducing V re v The complicated 

285 behavior discussed above convinces us that sometimes numerical results cannot necessarily help us if 

286 we do not know the whole truth from analytical results. 



Figure 4. Graphs ofC\(X) when T>\ is fixed, but we increase T> 2 - 


287 2.3. A comparison with Goldman-Hodgkin-Katz equation for VW 


In this section, we first recall the GHK equation ([20,22]), which relates the reversal potential with 
the boundary concentrations and the permeabilities of the membrane to the ions. If the membrane 
is permeable to only one ion, then that ion's Nernst potential is the reversal potential at which the 
electrical and chemical driving forces balance. The GHK equation is a generalization of the Nernst 
equation in which the membrane is permeable to more than just one ion. The derivation of GHK 
equation assumes that the electric field across the lipid membrane is constant (or equivalently, the 
electric potential (p(x) is linear in x in the PNP model). Under the assumption, the I-V (current-voltage) 
relation is given by 


I = V t 4 D k 


k=l 


n ~ he ZkV 

1 - e z * v 


For the case where n — 2 and z-\ — —Zi — 1, the GHK equation for the reversal potential is 


V GHK 
v re v 


( P ) = l n 


r + pi 
l + pr' 


( 8 ) 


The assumption that the electric potential c p(x ) is linear is not correct when applied to channels in 
proteins. That is because proteins have specialized structure and spatial distributions of permanent 
charge (acid and base side chains) and polarization (polar and nonpolar side chains). Experimental 
manipulations of the structure of channel proteins show that these properties control the biological 
function of the channel. The GHK equation does not contain variables to describe any of these 
properties and so cannot account for the biological functions they control. A linear <p(x) is widely 
believed to make sense without channel structure presumably, in particular, where Qo = 0. However, 
this is not correct, either. It follows from (7) for Qo = 0 that the zeroth order in e approximation of the 
reversal potential in this case is, 

V rev (0,p) = ^\n 1 -. (9) 

p +1 r 

288 Figure 5 compares V rev (0,p) in (9) with from the GHK-equation in (8). It can be seen, from the 

289 left panel that when / and r are close (for example L = Cq/ = 20 mM, R — Cor — 50 mM), then the two 

290 curves have almost the same behavior. However, when we reduce L from 20 mM to 4 mM, then the 

291 right panel shows a significant difference between the two graphs. 
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Figure 5. V rel ,(Q = 0,p) vs Vjfy K (p): The left panel for L = 20 mM and R = 50 mM; the right panel for 
L = 1/i/M and R = 50 mM. 

In Figure 6 , we arrange a simple numerical result for the case where Q 7 ^ 0 to compare the graphs 
of V r ev(Q,p), obtained from (7), for various values of permanent charge Q. We consider L = 20mM, 
R — 50 mM, and 0 < p < 5 for some values of Q, i.e., Q — 0 M, 1M, 10 M. 



Figure 6. V rev (Q,p) with various values of permanent charges. 


2.4. Profiles of relevant physical quantities. 

It follows from (1) and (2) that for any given Q, once a solution (A, V) are determined, all the other 
unknowns can be determined. We consider the dimensional form of quantities, and fix ( Q, L, R, 'D \, V 2 ) 
to numerically investigate the behavior of Cp(X) and <t>( X) throughout the channel. Figures 7 and 
11 graph the cases with small permanent charge Q — 0.1 mM when L — 20 mM, R — 50 mM, T>\ = 
1.334 x 10~ 9 m 2 /s, and £b = 2.032 x 10~ 9 m 2 /s. In this case, we obtain J — — 72.7387Mmol/(s N^) 
and Vrev — —4.4820 mV. 
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Figure 7. The functions Q-(X) (left) and < J>(X) (right) until Q = 0.1 mM. 

Furthermore, Figures 12 and 13 show graphs of concentrations, electrical potential, and 
electrochemical potentials versus X, where L = 20mM, R — 50mM, Q — 2M, and diffusion 
constants are the same as previous one. In this case, we obtain J — — 11.3146Mmol/(s N,i) and 
V rev = -19.5527 mV. 

Remark 3. We end this section with a few of remarks on some important features captured in the 
figures 7, 11, 12, and 13. It follows from Nernst-Planck equation that }i' k (x) has the same sign as that of 
f(jt(l) — /ijt(O) or the opposite sign as that of f t; in particular, }ik(x) is always monotonically increasing 
or decreasing. For zero-current situation, the reversal potential depends on ALL other parameters; 
and so it seems that it would be hard to make general conclusions about }i k {x), for example, about its 
monotonicity. This is not true. In fact, in Section 2.1, we have concluded that the sign of zero-current 
flux / is the same as that of L — R, and hence, }i’ k (x) has the opposite sign as that of L — R. For the 
cases considered in this part, L — 20 mM < R — 50 mM, one has / < 0, independent of the value of 
Q. Therefore, }i'Ax) > 0 for k — 1,2, and hence, }i k (x)’s are increasing for zero-current situation when 
L < R, independent of Q, as shown in Figure 11 for Q — 0.1 mM and in Figure 13 for Q — 1 mM. 
On the other hand, as one changes the value of Q, the profiles of c k (x)’s and <p(x) may vary from 
monotone to non-monotone, as shown between Figure 7 for Q = 0.1 mM and Figure 12 for Q = 1M. 
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Figure 8. The functions jq(X) and ;< 2 (X) are increasing for Q = O.lmM. 
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Figure 9. The functions Ci(X) and C 2 (X) (left) and the function 4>(X) (right) for Q = 1M. 



Figure 10. The functions }i\{X) and ft 2 (X) are increasing for Q = 1M. 

where L — 20 mM, R — 50 mM, Q = 2 M, and diffusion constants are the same as previous one. 
In this case, we obtain J = — 11.3146Mmol/ (s N/\ ) and V rev — —19.5527 mV. 

Remark 4. We end this section with a few of remarks on some important features captured in the 
figures 7, 11, 12, and 13. It follows from Nernst-Planck equation that }i[(x) has the same sign as that of 
^j-(l) — ^jt(O) or the opposite sign as that of /j.; in particular, }ik{x) is always monotonically increasing 
or decreasing. For zero-current situation, the reversal potential depends on ALL other parameters; 
and so it seems that it would be hard to make general conclusions about fi^ifx), for example, about its 
monotonicity. This is not true. In fact, in Section 2.1, we have concluded that the sign of zero-current 
flux / is the same as that of L — R, and hence, has the opposite sign as that of L — R. For the 
cases considered in this part, L — 20 mM < R — 50mM, one has / < 0, independent of the value of 
Q. Therefore, > 0 for k = 1,2, and hence, y/ c (x)'s are increasing for zero-current situation when 
L < R, independent of Q, as shown in Figure 11 for Q — 0.1 mM and in Figure 13 for Q — 1 mM. 
On the other hand, as one changes the value of Q, the profiles of cy (x)'s and fix) may vary from 
monotone to non-monotone, as shown between Figure 7 for Q = 0.1 mM and Figure 12 for Q = 1M. 
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Length X in nm 

Figure 11. The functions /<i(X) and J< 2 (X) are increasing for Q = 0.1 mM. 




Figure 12. The functions Cj(X) and C 2 (X) (left) and the function <I>(X) (right) for Q = 1M. 



Figure 13. The functions /<i(X) and f< 2 (X) are increasing for Q = 1M. 

332 3. Current-voltage and current-permanent charge behaviors. 

333 Ionic movements across membranes lead to the generation of electrical currents. The current 

334 carried by ions can be examined through current-voltage relation or I-V curve. In such a case, voltage 

335 refers to the voltage across a membrane potential, and current is the flow of ions through channels 

336 in the membrane. Another important data is current-permanent charge (I-Q) relation. Dependence of 

337 current on membrane potentials and permanent charge is investigated in this section for arbitrary 

338 values of diffusion constants. 
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To derive the I-V and I-Q relations, we rely on [13] where the authors showed that the set 
of nonlinear algebraic equations is equivalent to one nonlinear equation for A, defined in (3). All 
other state variables such as fluxes, profiles of electric potential <p(x) and concentrations c k (x) can be 
obtained from the special variable A. It is crucial to realize that this is a specific result, not possible 
in many cases. One can only imagine that the resulting simplification produces controllable and 
robust behavior that proved useful as evolution designed and refined protein channels. The reduction 
allowed by this composite variable can be postulated to be a "design principle" of channel construction, 
in technological (engineering) language, or an evolutionary adaptation, in biological language. In 
particular, the current I can be explicitly expressed in terms of boundary conditions, permanent charge, 
diffusion constants, and transmembrane potential in the special case that allows the determination of 
A. In what follows, we derive flux and current equations - when diffusion constants are involved as 
well - in terms of boundary concentrations, membrane potential, and permanent charge. The I-V, I-Q, 
J-V, and J-Q relations are investigated afterward in Section 3.2. 


3.1. Reduced flux and current equations. 

It was shown in [13] that the singular orbits of BVP (9) can be reduced to the algebraic equation 


rj In ^—— — N — 0, 
S a -tJ 

where B — l — A + r, and S a , S b and N were defined in (4), and. 


V = Qo - 


Qo 

In 

111 Ar 


V + ln 


l(S b ~Qo) \ N 
r(S a -Q 0 )J Inf 


Once A is solved from (1), we can obtain the flux densities and current equations as follows, 
Jk := J fc (V,/,r,D 1 ,D 2 ) =3D k (l-A)(l + (-l) fc ^-), fork = 1,2, 


I := 


I(V,l, r, D lf D 2 ) =Ji ~ h = 3(Z - A) (o a - D 2 — (D x + D 2 )). 


( 1 ) 


( 2 ) 


(3) 


For any given (l, r, D\, D 2 , Qo, V) there exists a solution for the flux / and current I. The numerical 
results in the next section give us more information on "current-voltage" and "current-permanent 
charge" relations. 


3.2. Current-voltage and current-permanent charge relations. 


Dependence of current on diffusion constants. Now we reveal a feature of the theoretical results that 
is not intuitive. Suppose that (l, r, Q) is given (V is still free and can have any value!). It follows from 
(4) for the definition of N that there exists an A so that N — 0. It consequently follows from (2) that, if 


V = In 


5(S«-Qo) 

A(S b - Qo) 


then // 


0. Therefore, from (3), I 


3(1 — A)(D i — D 2 ), which implies. 


For some fixed parameters (l, r, V , Q), the sign of I depends on the sign ofD\ — D 2 . 


I-V curves and I-Q curves. Figure 14 is a numerical simulation from (1) and (2) of the I-V curves for 
several values of Q with T>\ — 1.334 x 10 -9 wz 2 /s and X> 2 — 2.032 x 10~ 9 m 2 /s. One may suspect, based 
on the numerical observations, that the value of current X, obtained from (3), is unique for any V, and 
is monotonically increasing in V. However, this may not be correct, in general. This is important since 
the opening and closing properties of channels might be thought to arise from non-unique solutions 
([9,10]). On the other hand, our numerical experiments shows that: 

I-Q curves are not monotone in general. 
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Figure 14. The function I = X(V) for L = 20 mM and R = 50 mM. 

Note that the equation (6), in dimensional form, gives, 

vJx)cm dX=M0) - Mt) ' t = hl 

The sign of is determined by the boundary conditions, independent of the permanent charge Q. 
Nevertheless, as expected and seen in Figure 15, the magnitudes of Jp s, and consequently, the sign and 
the size of the current X do depend on the permanent charge Q in general. The numerical observations 
in Figure 15 indicate that 

(i) there exists some V* such that, for V > V*, I ma x{Qo ) has a unique maximum; 

(ii) there exists some V* such that, for V < V*, X m ,„(Qo) has a unique minimum. 

In particular, I-Q curves are not monotone in general. 




Q in M Q in M 


Figure 15. The function X — X(Q) with T>\ < XT-' The left panel for L = 20 mM and R = 30 mM; the right 
panel for L = 30 mM and R = 20mM. 


Besides, we claim based on numerical observations (not proven though) that there exists 
V(X>i,X> 2 ) — min{V*, |V*|}, such that 

(i) for any given V where | V | > V, the corresponding current X is non-monotonic in Q, but 

(ii) for any V where | V| < V, the corresponding current X is monotonic in Q. 

In particular, it can be seen, in Section 3.3, that current is monotonic in Q for V = 0. Tn the end, we 
would like to mention that the diffusion constants affect the values V* and V* above. 
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3.3. Zero-voltage current. 

The different permeability of the membrane determines the zero membrane potential (voltage) to 
different types of ions, as well as the concentrations of the ions, the permanent charge, and the shape 
of the channel. Denote current I(V; Q), and the fluxes JifV; Q), for k — 1,2, to include the dependence 
on Q too. We call I (0; Q) the zero-potential current and ] k (0; Q) the zero-potential fluxes respectively when 
V = 0. For any given value of membrane potential V , approximation formulas for the current I(V;Q), 
for small and large values of permanent charge Q, are provided in [29] and [48] respectively. 

For simplicity, we consider the case where a — 1/3 and fl = 2/3. It follows from [29] that for 
small values of Q, applying V = 0, zero-potential current I s ( 0; Q), and zero-potential fluxes J|(0; Q) 
(in dimensionless forms as mentioned in Remark 2) are. 


I s (0;Q) =(/ 

7k(0;Q)=(/ 


- r) (D 1 - Do) - + 0(Q 2 ), 


-r)D k + (-!)* 


2(2Z + r)(/ + 2r)ln' 
3 (l-r) 2 D k 


2(21 + r)(l + 2r) In - 


jQ + o(Q l ), k 


— 1/2. 


(4) 


For large positive values of Q — 2Qq, with v — 1 / Qq (where v is small), it follows from [48] that 
zero-potential current I 1 (v) — I 1 (0;Q) and zero-potential fluxes j[(v) — j((0;Q) are. 


where 




3 -d 2 / ! +r 

2 (vt+v^y 2 

l + r 


f(l,r)(l — r)v + 0(v 2 ), 




l + r 


f 2 (v) =6D 2 Vlr ^l ^ - |d 2 ‘ ^ r) (l - r)v + 0(v 2 ), 

Vl + V r 1 (VI + V r ) 


f(lr)=- 


2 Ir 


, 1 In / — In r , , 2 

+ ? + r ~2 /_ r ( l + r )■ 


(5) 


( 6 ) 


(VI + Vr) 2 

...It can.be readily.seen from equations (4) that, for small, values of Q, the zero-potential current 
/ (0; Q) is increasing m Q when f < r and is decreasing in Q if / > r. r 

But, for large values of Q, the zero-potential current 1 1 (0; Q) depends on Q in a much richer way. 

To state the results, we need the following lemma. 


Lemma 1. There are t\ and t 2 with 0 < t\ < 1 < t 2 so that f(l, r) > 0 for l/r G (ti,t 2 ) an d f(l,r) < 0 for 
l/r & [fi/ 121 - 


Note that 


±i\o) = + IM.) n _ r) , 

dv 2 KVl + Vr^ f + t/ 


It follows from (5) and Lemma 1 that, for large values of Q (small values of i'). 


(i) if l/r £ [t i, t 2 \ (so that f(l,r) > 0), then, for arbitrary Dj and D 2 , the zero-potential current I 1 (v) 
is decreasing in v (increasing in Q) when l/r € [t\, 1), and is increasing in v (decreasing in Q) 
when l/r G (l,f 2 ]; 

(ii) if l/r ^ [fi, f 2 ] (so that f(l,r) < 0), then. 


(a) for ^ + ‘ 7 +p > 0, the zero-potential current t 1 (v) is decreasing in v (increasing in Q) for 


l < r, and is increasing in v (decreasing in Q) for l > r; 
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(b) for ^ < 0, the zero-potential current 1 1 (v) is increasing in v (decreasing in Q) for 

1 < r, and is decreasing in v (increasing in Q) for l > r. 

Figure 16 illustrates some of the above conclusions. Besides, it suggests that the monotonicity 
of 1(0 ) holds for all values of permanent charge, not only for small or large values. We emphasize 
that the monotonicity of current I with respect to permanent charge Q is just true for zero membrane 
potential, i.e., V = 0. Indeed, one should recall from Section 3.2 that when V f= 0, then the current I is 
not monotonic in Q. 




Figure 16. The function T = T(Q) for V = 0: The left panel for L = 20 mM and R = 30 mM; the right panel 
for L = 30 mM and R = 20 mM. 

4. Conclusion. 

In this paper, we first recall the analytical results in [38] for arbitrary diffusion constants. To 
investigate the reversal potential problem for which the current is zero, we do numerical investigations 
based on the analytical results in [38], where many cases are studied analytically. We derive several 
remarkable properties of biological significance, from the analysis of these governing equations that 
hardly seem intuitive. 

Biophysicists are also interested in the relation of current-voltage (I-V), and current-permanent 
charge (I-Q), as well as reversal potentials problems. To do that, we first recall the analytical results 
in [13], for arbitrary diffusion constants, to drive the flux densities and current equations explicitly. 
One way to characterize channels is the current at zero electric potential, that is when V — 0, that 
has practical advantages. Since it is usually easier to measure a large current than a vanishing one, 
we analyzed this case, as well. Furthermore, we briefly study the special cases of small and large 
permanent charge for zero voltage case, based on the analytical results of [29] and [48] respectively. 
To bridge between small and large values of permanent charges, we numerically study I-V and I-Q 
relations for this case as well. 
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Physical Sciences-Collaboration Grants for Mathematicians #581822. 
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